Elastic properties of graphene flakes: boundary effects and lattice vibrations 
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We present a calculation of the free energy, the surface free energy and the elastic constants 
("Lame parameters" i.e, Poisson ratio. Young's modulus) of graphene flakes on the level of the 
density functional theory employing different standard functionals. We observe that the Lame 
parameters in small flakes can differ from the bulk values by 30% for hydrogenated zig-zag edges. 
The change results from the edge of the flake that compresses the interior. When including the 
vibrational zero point motion, we detect a decrease in the bending rigidity, k, by ^ 26%. This 
correction is depending on the flake size, A'^, because the vibrational frequencies flow with growing 
A'' due to the release of the edge induced compression. We calculate Griineisen parameters and find 
good agreement with previous authors. 
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I. INTRODUCTION 

Since its fabrication has become technologically 
feasible^, graphene has been in the focus of frontier 
research^Tii. One of its most celebrated properties are 
its massless low energy excitations^iS, ("Dirac fermions"), 
which emanate from the symmetries of the honeycomb 
lattice. The electronic properties of graphene flakes are 
quite different from bulk graphene due to the finite size 
and the presence of edgesii^. In particular, calculations 
suggest that the zigzag edges of graphene nano-ribbons 
(quasi Id) have two flat bands at the Fermi energy^'^° 
that introduce magnetisn>iir— . Recent theoretical stud- 
ies on zigzag edged graphene fiakes also confirm a ten- 
dency towards edge mag netismil"^. This gives addi- 
tional motivation for fabricating graphene based nano 
structures^"— for further studies. The fabrication of 
such structures with well defined edges still poses a con- 
siderable technological challenge. Therefore, only very 
few experiments with structures exhibiting zigzag edges 
have been reported22r— ; a detailed investigation of the 
edge physics still needs to be done. 

An increased interest in the elastic properties of 
graphene has developed recentl y^^'^"" — . This is, for in- 
stance, because experiments suggest that graphene sam- 
ples exhibit a corrugated structure^^i^"— ("ripples") 
even at relatively low temperatures. Their origin is 
thought to be due to residual elastic strain produced by 
the experimental preparation technique^. 

Also for elastic properties, edge effects can be highly 
relevant. For this reason, flake elastic properties are cer- 
tainly interesting in their own right. Namely, there is an 
intimate relation between the electronic structure and 
the atomic geometry of graphene. For example, the elec- 
tronic spectrum of a certain class of arm-chair graphene 
nano-ribbons is reported to acquire a spectral gap due to 
an edge induced lattice dimerization along the transport 




FIG. 1: Geometry of the NxN hydrogenated graphene 
fiakes(here A'^=4), that we have used for the density func- 
tional calculations. 



directiorti^,. 

In our study we investigate the free energy, the elastic 
properties, and the phonon spectrum oi N x N graphene 
sheets ( "flakes" ) as displayed in Fig. [1] using the den- 
sity functional theory (DFT). We show, that the different 
chemical nature of C-C bonds at the hydrogenated edge 
as compared to the bulk leads to an nearly homogenous 
compression, i.e. strain. As a consequence, the average 
C-C-distance in a 3 x 3 fiake is reduced by a substantial 
amount, 0.3%; for comparison, strain as achieved in typ- 
ical pressure experiments does not usually exceed values 

^ 34.48. 49 

The presence of the surface induced strain leaves var- 
ious traces in the fiakes' interior observables. (a) The 
fiakes elastic constants, i.e. the Lame parameters, are en- 
hanced as compared to the bulk case. For isotropic strain 
in smallest fiakes (3 x 3), the (inverse) compressibility 
/i -I- A (precise definition see below, Eq. ([2])) increases by 
30%; for shear forces the increase is even bigger, almost 
a factor of 2. (b) Under bond compression the inter- 



2 



atomic forces typically increase, so that even the short 
wavelength vibrations, in particular the optical phonons, 
exhibit a "blue shift" of their frequencies with decreas- 
ing flake sizes. This flow can be seen in the variation of 
the Raman spectra with strain— i'^^'^"" — and can be de- 
scribed in the standard manner by Griineisen parameters. 
The values that we find here of Griineisen parameters, 
agree reasonably well with previous reports^i^i^. 

Even though one might suspect, that our topics have 
already been dealt with extensively in the literature^"—, 
a detailed investigation of the elastic properties of nano- 
fiakes is yet to be done; this refers, in particular, to an 
analysis of edge and finite size effects of hydrogenated 
zig-zag flakes which we perform in this work. 



II. FREE ENERGY OF HOMOGENOUS, 
PLANAR FLAKES 



We consider the free energy of a graphene flake as de- 
picted in Fig. [T] with a homogenous C-C-distance, d, as 
the sum of all its bond energies: 



N, 



(1) 



where A'i denotes the number of internal C-C-bonds with 
an associated binding energy VP, N^. denotes the num- 
ber of edge located C-C-bonds with energy -0 and -0'^ 
includes the corner contributions, where Nc is the num- 
ber of bonds linking the corner atoms, -/Vc=6 in Fig. [1] 
The binding energy per C-C-CH edge group is close to 
2ip but not identical to it. For instance, tp also includes 
corrections of internal bonds, that still "feel" the pres- 
ence of the surface. Similarly, Nc^'^ is approximating the 
binding energy of the corner groups (two HC-CH groups 
and two C-CH-C groups). 

We mention, that the representation ^ is slightly sim- 
plified in the following sense. In general, the boundary 
(shape) of a given flake, e.g as depicted in Fig. [U does 
not share the hexagonal symmetry of the honeycomb lat- 
tice. For this reason, in flakes with a fully relaxed atomic 
structure bond lengths and bond angles are not strictly 
all the same. Our DFT-calculations indicate, however, 
that such distortions, though clearly detectable, give only 
small corrections to those phenomenological parameters 
that we are mostly interested in. 

The continuum theory of 2c?-membranes has been de- 
vised for an inhomogeneous flake with neighboring bonds 
exhibiting slowly (in space ) varying bond distances, d{r), 
and angles. In this formulation the elastic energy is rep- 
resented by the functional^, 

E ^ ^ I d^iW^hf + f d^r (u^^+uyy)' 

^2!/'' + • (2) 

The flake coordinates are given with respect to a planar 
reference state with area A, that lives in the r=(a;,?/)- 



plane; accordingly, the in plane coordinates constitute 
the displacement vector, u(a;, y), that measures the trans- 
lation of each membrane point (x, y, z) with respect to 
the reference state. The out of plane distortions deflne 
the height field h{x,y); for the planar case h—0. u and 
h together constitute the strain tensor (i,j = x,y), 
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It is clear that the form of Uij depicted in Eq. ^ is sym- 
metric construction, with respect to the spatial deriva- 
tives of u ensure its invariance under in plane rotations 
by 90° which correspond to u cx {y, — x). 

The total elastic energy is a sum over contributions 
which resemble local oscillators in the membrane plane. 
The first term is proportionate to the curvature V^/i and 
introduces the bending rigidity k. It describes the en- 
ergy cost for bending the membrane without changing 
the bond lengths or in plane bond anglesi^^ The Lame 
parameters, A and /i, appearing in the second and third 
term of Eq. ([2]) describe the in plane rigidity. 

For homogenous, planar membranes the elastic theory 
(p|) may be considered as a continuum approximation to 
([1]) which does not make explicit reference to boundary 
terms. Edges are accounted for only in the boundary 
conditions and (possibly) in a dependency of the Lame 
parameters on the position with respect to the edge. Usu- 
ally not included in ([2]) is the fact that this spatial de- 
pendency supports long range terms, ~ 1/flakesize. They 
modify the Lame parameters appearing in ^ even inside 
the flake's interior. 



Phenomenological parameters 

Isotropic strain: In order to illustrate the coopera- 
tive effect between surface and bulk, we consider an ex- 
pansion of ([T]) in terms of the variable e={d — do) /do; e 
quantifies the strain inside the flake. The bulk free energy 
per bond has an expansion, 
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l^3e' + ^^,e\.. (4) 



where the bulk bond length do is to be determined at 
iVi.o ^ 00. The surface free energy may also be ex- 
panded about a minimum bond length, dg , but in general 
do ^ do- After all, in the limit A^i c — >■ 00, just the first 
term in Eq. ([1} contributes to the free energy per area 
and therefore do needs to minimizes 5*, only. Hence, we 
introduce the relative deviation of surface and bulk op- 
timal bond lengths, d—^do — do) /do, so that we have an 
expansion 



= 00 + ^02 (£+0)' + ^03 (£+0)' + ^04 (£+0)^ 
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where the coefficients in the second hne are defined in 
terms of the expansion the hne before. The elastic prop- 
erties of the flake are determined by the expansion pa- 
rameters ^2,3,4, V'l, 2,3- 

At any finite value of Ni^c, optimization must also in- 
clude the boundary (i.e. surface) terms and therefore 
the optimal value of e, ejv, is non-vanishing in this case; 
specifically, 



Eat 



*2 Ni 
*2A^i' 
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In order to calculate the feedback of this shift into the 
elastic parameters, we expand F in the vicinity of its 
minimum, ejvi to the fourth order in e. Recalling that 
this corresponds to a strain u(x)=ex we can compare the 
result with Eq. ^ and thus find: 
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(7) 



The first term in the rhs-expression ([T]) simply accounts 
for the separate, additive contributions of bulk and sur- 
face (i.e. edge) free energies. The edge contribution, that 
appears here, could formally be accounted for in a gener- 
alized version of Eq. ([2]) where one adds a boundary term. 
Similarly, by allowing for a dependency of the Lame pa- 
rameters on strain itself, one could also include additive 
an-harmonic effects, second bracket first two terms. In 
either case, the phenomenological parameters are univer- 
sal in the sense, that they are the same for each flake size 
and geometry. 

The interesting pieces are the terms in sn, which are 
mixing surface and bulk parameters: vPsEat ~ ^'1^3/^2- 
They encode the "cooperative" effect between boundary 
induced strain and bulk anharmonicities. It is due to 
them, that the flakes elastic parameters need to be ad- 
justed in principle for every geometry separately. 

Shear strain: An analogous analysis as for the 
isotropic strain also applies to shear forces. The expan- 
sion is even in the shear strain u = (0, Egx): 



(8) 
(9) 



The new expansion parameters, ipiji — 2,4,.. ., are, 
in general, dependent on the flake geometry. Again, by 
comparing to Eq. Q we flnd 



(10) 



^ = -f + Y^e; ( *4 + 



Here, surface and bulk free energies give strictly additive 
contributions, and a cooperative effect does not emerge. 



III. DENSITY FUNCTIONAL CALCULATIONS 
A. Method 

In Fig. [T] we display the geometry of the N x A^- 
graphene flake that is employed in our calculations: Ni — 
(iV-l)(3A^-l), iVo = 8(iV-l). Electronic structure cal- 
culations have been performed for a given atomic config- 
uration (C-C distance, flake geometry etc.) on the basis 
of the density functional theory as implemented in the 
quantum chemistry package TURBOMOLIi^. We are 
comparing GGA functionals (BPSe^SiSi^ PBE^SiS) with 
a hybrid functional (B3LYP~) and use a minimal basis 
set (SVP^-). Specifically, we are working at zero tem- 
perature and approximate the free (i. e. ground state) 
energy, Eq. ([T]), by the DFT estimate for the total bind- 
ing energy of the fiake: 



with. 



Fei(7V, d) := E,i{N, d) ~ Ei,,,iN). (11) 



Ef,,,{N) = NnE^ + NcEc (12) 



where i?H/c denote the DFT energies of a free charge 
neutral hydrogen/carbon atom and -/Vr/c denotes the 
number of hydrogen/carbon atoms in the flake. 



B. Results and Discussion 

Isotropic strain: A sequence of DFT calculations 
has been performed for A^=3 ... 9 and different values of 
the C-C distance, d. For each distance, Fc.\{N,d) has 
been calculated. In order to extract the expansion coeffi- 
cients of Eq. (P), ^{d),ip{d),^''{d),(j)%d), (j)'{d), we have 
performed five parameter fits on sets of raw DFT data. 
These fits were applied to three data sets consisting of 
{iV}=3, ... 7, {N}=3 ... 8 and {A^}=3 ... 9. The results 
for the surface, bulk and corner free energy have been dis- 
played in Fig. [21 The scatter between the fitting param- 
eters belonging to different data sets is relatively small, 
which illustrates the stability of the fit. 

The lattice constant of bulk graphene is estimated from 
the minimum position of ^'(d) Fig. [5J upper panel as 
(io=2.706ao, where ao=0.529A denotes the Bohr radius. 
Comparing this position to the minimum of the edge (sur- 
face) free energy. Fig. [2l center panel, dg=2.694ao, we 
find i)=0.44%. This indicates clearly the compression of 
the C-C bond length near the edge. The shift of the 
minimum position to lower values becomes even more 
pronounced near the corners, i.e. in ^'^{d), see Fig. [2l 
lower panel. 

Furthermore, we can perform a certain consistency 
check by evaluating the bond energies. We have a bind- 
ing energy ^(do)= — 5.22eV for bulk carbon atoms. The 
binding energy of _ff-atoms near edges (C-C-CH group) 
is approximately A£'|^«2-0((io)— 25'((io)~ ~ 4.6eV; when 
going to the corners (HC-CH and C-CH-C groups) we 
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FIG. 2: Bulk (^), surface (edge, ^) and corner {ip'^) free 
energy per carbon bond in the grapliene flakes, Fig. [TJ calcu- 
lated with density functional theory (BP86- functional). Data 
based on the evaluation of three sets of flake sizes ranging 
from {A''}=3 ... 9, and a 5 parameter fit to Eqs. (|Illip per d 
values. (ao:= 0.529 A.) 





-*o [eV] 


i*2 [eV] 


-i*3 [eV] 


BP86 




4g30^±O.082 


128.301±^-^i'^ 


B3LYP 


5.008±«*^ 


47.723*°-^'^^ 


154.59±^i-^^« 


PBE 


5373±o.oo4 


45.57±°-« 


187.309±2*-*** 



TABLE II: Bulk free energy coefficients as defined in Eq. 
These coefficients are extracted from fitting Eq. Q to the 
data in Fig. [21 upper panel. 





-^0 [eV] 


V-i [eV] 


iV2 [eV] 


-iV3 [eV] 


BP86 


75^5±o.oo6 


285±o.oi5 


4493±o.43 


112.745=^^*-^^* 


B3LYP 


7 i87±''.oo6 




45.586±2-^^'* 




PBE 


7g3±o.oi 


0.308±°«" 


47.503*^'"*^ 





TABLE III: Edge (surface) coefficients as defined in Eq. ((5|. 
These coefficients are extracted from fitting Eq. ((SJ to the 
data in Fig. (2] middle panel. 



the slope and perhaps also the sign of the tvifo functions 
has converged, already. Under this assumption we may 
conclude, that both amplitudes flow closer to zero val- 
ues when e increases. This behavior is compatible with 
the simple expectation, that the main effect incorporated 
in the l/iVj^e-corrections is the discreteness of the flake's 
electronic spectrum with level spacings Aj c for bulk and 
surface modes. With increasing e the bandwidth de- 
creases and so do Ai e and 0'"°. 





do [A] 


dS[A] 


BP86 


^ 432±o.oo2 


1 427*"''-'''-'- 


B3LYP 


1.426±»*i 


I.42I±"-*2 


PBE 


^43^±o.ooi 


^ 426±o.oo5 



TABLE I: Minimum C-C bond length as extracted from 
bulk free energy and surface free energy correspondingly (see 
Fig. [21 upper and middle panel). The distance in do and d^ 
leads to a nearly homogenous pressure on the flake that modi- 
fies elastic and electronic-structure properties. Data is shown 
for three different functional used in DFT calculation. 



have AE^Kiil;'^{do)—'^{do) « —5.3 which is roughly con- 
sistent with Ai^lj, as it should be. 

To obtain also the other phenomenological parameters, 
a second (polynomial) fit of the traces ^(d), ip{d), Fig.|5J 
according to Eqs. (|4l5p has been performed; all fitting 
parameters are summarized in Tab. Illlll and Hill 

When fitting the raw data to get 4', V', V''^ the terms in 
l/iVi.c could not be neglected for the system sizes, that we 
considered. The corresponding amplitudes are displayed 
in Fig. 131 Unlike it was the case with the previous data. 
Fig- m the amplitudes (/)''° of the l/A^i.o-corrections still 
exhibit a considerable variation with increasing system 
size, which is probably due to even higher order terms 
that have been neglected in the expansion Eq. In- 
terestingly, while the magnitude of 4>^'°{e) is still shifting 



Shear strain: A largely analogous method as was 
adopted for the isotropic strain, has also been applied 
for shear forces. In this case, the convergence of the 
DFT calculations turned out to be considerably more 
difficult, so that the investigated system sizes range from 
{A''}=3 ... 7, only. From our fitting procedure we could 
determine the response of the bulk and surface energy to 
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FIG. 3: Dependency of the amplitudes i^''"' describing the 
corrections in 1/Ne,i to the binding energy l\Fs\{N). Data 
were obtained by the 5 parameters fit, already underlying the 
traces shown in Fig. [2l 
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FIG. 4: Change in the free energy when applying a shear 
strain in graphene flakes. Upper panel shows the change in 
bulk free energy per interior C-C bond; lower panel exhibits 
change in surface (edge) free energy per C-C bond. Data sets 
are for system sizes {N}=3. . .7 and extracted from an expres- 
sion analogous for shear to Eqs. . The lines indicate the 
polynomial fit according to Eqs. (|8l9p with parameter given 
in Tab. [Wl 



[eV] 


i*2 


24 * 4 






BP86 








5237±60% 



TABLE IV: Bulk and surface shear free energy coefficients as 
extracted from fitting Eq. ((8)1 (Eq. (|9])) to the data in Fig. |4] 
upper panel (lower panel). 

the shear strain, eg, as shown in Fig. S) The parameters 
entering Eqs. (I8I9P can be extracted and are listed in 
Tab.llVl 

IV. FLAKE ELASTIC PROPERTIES 

In the previous subsection, the focus was on the be- 
havior of the free energy on the flake size under isotropic 
and shear strain. In this section, we discuss and illustrate 
what our previous findings imply for the elastic proper- 
ties of a single flake with a fixed size, iV. Partly, we are 
considering the same set of data again, but now plotting 
observables directly for iV fixed. 

1. Homogenous isotropic strain 

Fig.[5]shows, how the excess energy per unit cell grows 
under increasing strain for different flake sizes N. It is 
readily seen from this plot, that there is a shift of the 
equilibrium lattice constant d{N) to smaller values. In 
the light of the previous section, this shift is the expected 
consequence of the surface induced strain en- The inset 
shows, the scaling with Ne/Ni. 

In addition, we also extract the flake elastic constant 
li+X from the parabolic shape of the curves, Fig. [5] To 




FIG. 5: Excess energy AE per unit cell generated by rescal- 
ing all bond length, d, (homogeneous, isotropic strain). AE 
exhibits flow of the equilibrium bond length with the linear 
flake size A*'. Main panel: AE over bond length d employed 
in the simulation, d is measured relative to the bulk bond 
length. The lines serve as a guide to the eyes. Inset: Extrap- 
olating the equilibrium bond length d{N) into the bulk limit: 
do = 2.707 ±0.001. 



this end, we replot the data in Fig. [5] left, so as to high- 
light the curvature and its strain dependency. On the 
basis of Eq. ([7]) we can conclude, that the offset of the 
curves is a consequence of (a) the presence of the surface 
and the extra energy required for its compression (term 
Nc/Ni in Eq. ([7|)) and (b) the feedback of the surface 
strain sn into the bulk C-C distance. Extrapolating the 
zero strain values into the limit, iV^oo, we recover the 
previously derived bulk limit value. This check is dis- 
played in Fig. [SI right. The plot also reveals, that the 
deviation of elastic parameters from their bulk values in 
small flakes may not be very small. For our smallest 
flakes, N=3 it reaches almost 30%. 

Additional information can be extracted from Fig. [S] 
left, about anharmonicities which manifest themselves in 
the slope of the curves displayed. This pre-factor of the 
an-harmonic term (linear in e) in Eq. ([7]) admits the fol- 
lowing interpretation. The slope changes with increasing 
N since the contributions of the surface ('03-term) and 
the surface induced bulk compression (^'4ejv-term) di- 
minish. A non- vanishing value of ^'3/12 for the slope 
will remain however even in the bulk limit. 



2. Shear strain 

Following the same strategy as we did before with 
Fig. [6l we plot in Fig. [7] the excess energy AEg induced 
by pure shear strain, u(x) = es(0,a;). Again, the plot 
emphasizes the curvature in this quantity, /x, and how 
it evolves with the flake size. Since AEg is even in the 
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FIG. 6: Estimating the sum of the Lame parameters fi^i + 
Aci and the boundary correction (offset of traces) from AE 
displayed in the previous Fig. [5] Left: Data for curvature 
exhibit a slope which indicates the linear dependency of the 
Lame parameter Ad on strain. (Linear terms in /iei do not 
appear, see Fig. (9]) Right: extrapolating the curvature at 
£ = into the bulk limit. 





do [A] 


Aid + Aei[eV] 


Mcl[eV] 


f 


Y[N/m] 


BP86 


1 432±o.ooi 


70.715±°"" 




0.162 


356.23 


B3LYP 


^ 427±o.ooi 


71.21±oi2 








PBE 


^43^±0.002 










prcv. 

calc. 


1.42[37.66] 
1.41[43] 
1.45[68l 


66.5711371 


49.45[37| 


0.1491671, 
0.173[311 
0.16f68l 
0.311391 


3461371 
3071411 
3451671 
3361681 
312[391 


cxpt. 
(Graphcnc) 










342[34| 


cxpt. 
(Graphite) 


1.421[69] 
1.422[711 






0.1651701 


371[711° 



"assuming graphene thickness 0.335 nm. 

TABLE V: Comparison of C-C-bond distance in bulk 
graphene, elastic constants , Poisson ratio(!/) and Young's 
modulus(Y) as extracted from Fig. (|5I6I7|I respectively by ex- 
trapolating the values in bulk limit (TV — !> 00) with previous 
works. Data is shown for three different functionals we used 
in DFT calculation. 
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FIG. 7: Estimate for Lame parameter /iei determined from 
the excess energy AiJs per unit cell under pure shear strain 
with strength £s. (Procedure similar to previous Fig. [G]) Left 
panel: dependency of curvature of Ai5s(es) on the linear flake 
size A'^. Due to mirror symmetries of the unit cell, linear 
corrections do not appear for the shear parameter ^oi- Right: 
extrapolating the curvature into the bulk limit. 



We present results from an additional DFT study, 
where we investigate the transverse stiffness of the 
graphene flake that gives rise to the elastic parameter 
K. To this end we employ the following strategy. Each 
flake has a center pair or center ring of carbon atoms, see 
Figs. (|ll8p . To create a transverse probing held /i(r), we 
lift the center atoms by the distance /iq over the refer- 
ence plane. After this, the atomic structure of the flake is 
relaxed under the constraint, that the set of edge atoms 
(H-atoms and edge C-atoms) can move only within the 
reference plane; edge atoms cannot shift in /i-direction.— 
In this way, a flake is equipped with a single ripple while 
at the same time the associated strain fleld u(x) remains 
negligibly small. In order to estimate the integrated cur- 
vature we numerically compute the bi-variate function 
which interpolates the scattered data values(/i(r)-fleld) 
at any predeflned smooth mesh. We then use this inter- 
polated function to perform the second order numerical 
derivative at any arbitrary precision. 

Fig. displays how the excess energy /S.Eh associated 
with the ripple grows with the increased integrated cur- 
vature. 
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shear strain, only positive values of Es are given. Also, 
for the same reason an-harmonic terms exist only in the 
quartic order, so that the displayed data traces have zero 
slope. Similar to the previous case of isotropic strain, we 
also witness here a very strong dependency of the elas- 
tic constant on the flake size. In fact, for shear strain 
it reaches almost 70% for the small system sizes that we 



(13) 



The increase is linear, as expected from Eq. ([2]) with 
a slope that is only weakly dependent on the flake size, 
see inset Fig. [9l this implies that nonlinearities remain 
small as long as the ratio of the ripples amplitude and 
wave length , /iq/L, does not exceed ~ 5%. The bending 
rigidity thus found is n^i = 1.24 eV which is well consis- 
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FIG. 8: Buckling flakes of N=5, 6 with different central con- 
figurations of carbon atoms. The atomic configuration of C 
atoms is relaxed under the constraint that the center atoms re- 
main at a given height ho above the ground plane, while edge 
atoms (H and C) remain sitting within this plane {h = 0). 



tent with the value l.leV obtained by Fasolino, Los and 
Katsnelson^. 
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FIG. 9: Estimating the Lame parameter k determined from 
the excess energy AE of a bulging flake with maximum height 
at ho over the unperturbed (flat) plane (see Fig. |8|. Main 
panel: Change of energy AE with the integrated curvature 
I„ = f d?v (Ah(r)f for different flake sizes iV. Inset: The 
ratio K = AE/Ik depend on A'' due to the effect of edge 
compression. 



Notice, that there is a significant scattering of almost 
20% in earlier theoretical estimates for k and derived 
quantities, see Tab. 2 in Ref. Discrepancies ap- 



pear because different theoretical techniques are being 
employed, e.g. empirical potentials^ and density func- 
tional theory, but also because of modeling artefacts. For 
example, extracting k from the elastic energy of carbon 
nano-tubes (radius R) requires a very careful extrapola- 
tion in 1/R. If sub- leading terms are ignored, there is a 
pronoimccd tendency to overestimation, e.g. k = 1.46 eV 
in Ref. 67. These authors used nano-tubes with smaller 
tube radius, hence bigger curvature, where nonlinear ef- 
fects become important. We can check the bending rigid- 
ity usin g th e same curvature in Fig. [5] (inset) as reported 
in Ref. |67| and find a reasonable agreement with their 
value. 



V. ZERO POINT MOTION 

In this chapter we extend our analysis of flake elastic 
properties and take also the zero point motion of the 
atom cores into account, that constitute the hexagonal 
lattice. Now, the free energy acquires a second term, 

F(iV,d) =Fei(iV,d) + i^vib(Af,d) (14) 



with 



i^vib - (15) 



where p labels all the flake's vibrational modes. The 
vibrational excess energy associated with stretching the 
flake reads 



where uj{N,d{N)) denotes the vibration energies in the 
absence of strain and d{N) the equilibrium bond length, 
see inset Fig. [5l Also Ai^vib can be expanded in terms of 
the slow elastic modes, 



AFvib[/i,u] 



(17) 

with expansion parameters 7/1, 7„ that represent averages 
of Griineisen parameters over all vibrational modes. In 
a two-dimensional sample^ with a mirror symmetry one 
expects 7^^ — 7^^ — 0; the change in phonon frequen- 
cies should be even in the shear strain, Eg. Combining 
Eq. (jl7p with an expansion of Ai^oi in full analogy with 
Eq. (m and after completing the square, we flnd. 



J 



g1 vib 



{N,d) 



Kc\ + Ih 



d^v{V^hf 



d^Y 



In 



Mel + Acl 



^ / d^r(u,. 



"xy 



(18) 
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FIG. 10: Left: Dependency of AF^n, on isotropic strain 
for system size N—6. The slope defines 7„ as in Eq. (|17p . 
Right: Dependency of AF^n, on integrated curvature, 
/ d'r {Ah{r)f, as in Fig.H 



For clarity, we have indicated in this expression the bare 
electronic coefficients (i.e. with frozen atomic cores) by 
Kci, /.toi, Aoi. Likewise, the displacement field u(x) is de- 
fined with respect to the optimum fiake geometry ignor- 
ing vibrational terms. 

In this way we can observe two facts, (i) Vibrations 
modify the bare transverse stiffness KcI in Eq. (fTSi) . k = 
Kei + Jh- (ii) Vibrations also affect interatomic distances. 
The effect can be understood as an effective strain, which 
stretches the bare C-C distances: 7"/(^oi + Aci). (iii) The 
change in the C-C bond lengths eventually feeds back 
into all elastic coefficients. Therefore, in a more com- 
plete treatment of higher order terms also modifications 
in /Xoi, Aei would occur. 



A. DFT calculations of the Griineisen parameters 

In order to estimate the Griineisen parameters, 7/t, 7«, 
we should calculate the vibrational spectrum of flakes 
with and without applied strain. To this end we adopt 
the following procedure. For every flake size, N, we find 
the atomic geometry with the optimal electronic energy, 
see e.g. Fig. [T] This constitutes the set of freely relaxed 
"parent states" . The relaxation ensures the Hessian, that 
characterizes interatomic forces, to become a positive def- 
inite matrixjS 

In the present study we focus on the impact of phonons 
on the bulk elastic constants. There we may eliminate 
contributions of surface vibrations by assigning an infi- 
nite mass to the surface H and C atoms. Other than this, 
the calculation of vibrational modes and frequencies for 
the relaxed flake is a standard procedure^^i^. 

Thereafter, each parent state thus obtained is used in 
order to create two new families. The first family is con- 
structed to obtain 7„. It derives by changing the bond 
length of edge C-C-pairs by a factor of 1 -I- e keeping 
all atoms still inside the base plane (h—O). For each 
value £ the internal C-atoms are relaxed and the vibra- 
tional spectrum together with the average strain field, 
lui^) — J dr{uxx+Uyy), are recalculated. In this pro- 
cess it is important to have edge atoms at infinite mass. 





BP86, N=6 


prev. calc. 


expt. 




-0.055 


- 


- 


Ih [eV] 


-0.32 






7"/ (^el -1- Aol) 


-0.004 






7D [eV/A^] 


2.6 


2.7f66] 




7G [eV/A2] 


2.2 


2.0[75J 


1.99[48^ 



TABLE VI: Survey over the fitted Griineisen parameters ex- 
tracted from the data Fig. 1101 and Fig. [TT] respectively. 



This ensures that the fiake energy is in a (constrained) 
minimum so that all frequencies are real. 

In order to determine 7^ a second family has been con- 
structed. It consists of the buckled flakes. Fig. El that 
we have studied in the previous section in order to ex- 
tract Kei. Again, after assigning infinite mass to the edge 
atoms for each family member, the vibrational spectrum 
and the consecutive modification of the zero point energy 
can be calculated. 



B. Results and Discussion 

In Figs, [in] the change in the zero point energy, AFvib, 
is plotted over the integrated strain fields. The Griineisen 
parameters are given by the slope near zero strain; their 
numerical values are listed in Tab. I VII For a discussion 
of our results we first recall, that the vibrational spec- 
tral density of states of the carbon sheet has a strong 
peak in the optical frequency regime, c.f. Fig. Illl near 
1600 cm^^. It is the "G-peak", that corresponds to an 
in-plane mode, where neighboring atoms vibrate in op- 
posite direction as depicted in Fig. [T2] (right). This mode 
gives the dominating contribution to the total zero point 
energy. Another significant contribution comes from the 
frequency range 500— 1000 cm^^, where one observes the 
mixing of out-of-plane modes with in-plane modes. A 
third important mode is the "D-peak" near 1350 cm^^, 
that reflects the breathing mode as shown in Fig. [T^] 
(left). This mode is particularly interesting when study- 
ing finite size(edge) effects in graphene fiakes. The reason 
is that, in bulk graphene the coupling of the D-peak to 
the electromagnetic fields is suppressed, since the D^h- 
symmetry of the hexagonal unit cell remains intact and 
inhibits the formation of the dipole moment. By contrast, 
in flake with an overall symmetry that is lower than D^h, 
the D-peak is observable with a strength proportional to 
the inverse fiake size. 

Therefore, we can understand the sign of 7/1 as a con- 
sequence of a softening of these modes in the sample re- 
gions with non-zero curvature (A/i(r))^. Similarly, 7.^ is 
negative, indicating the increase of the atomic oscillator 
frequency that occurs when the interatomic distance is 
diminished. 

While the sign of jh.u was not unexpected, it is note- 
worthy that the vibrational contributions to the phe- 
nomenological material parameters are actually not so 
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FIG. 11: Left panel: Phonon density of states for the flake 
of size N=6. Data sets are shown for different C-C bond 
distances. Upper right panel shows the How of D-peak with in- 
plane strain. Lower right panel shows the variation of Raman 
G-peak with strain, where ojq is the Raman band frequency(~ 
2680 cm~^)^. The corresponding vibration modes are shown 
in Fig. [m 



small. The bare electronic bending rigidity, k^i is re- 
duced by as much as 26% down to k—Kci + 7^=0. 88eV. 
Similarly, when expressing the effect of vibrations on the 
atomic lattice as an effective strain pushing the atoms 
to larger distances, then this strain reaches values up to 
0.4%. 

Here, we also calculate the Griineisen parameters as- 
sociated with individual modes (see Fig. [T^ . The right 
panel in Fig. [11] shows the flow of the Raman frequen- 
cies (upper half D, lower G) with the applied strain. The 
frequency decreases linearly (as described in Eq. (IT7)) ) 
with decreasing compression due to anharmonicity of the 
interatomic potential. The slope essentially estimates 
the Griineisen parameters, 7g and 7d for the vibrations 
shown in Fig. [121 Our results for 7g, 7d are consistent 
with the earlier experiments and first principal calcula- 
tions (see Tab. [VTl) available in literature^^i^i^ii^. 



VI. CONCLUSION 



The elastic properties of edge hydrogenated graphene 
flakes have been investigated employing the density func- 
tional theory (DFT) . Our study emphasizes the interplay 
between the edge and bulk properties which are mediated 
via long range elastic forces. 

Speciflcally, we are able to disentangle bulk, surface 
and corner contributions to the free energy together with 
the leading higher order corrections. The binding energy 
per surface (edge) bond (7.5eV) is roughly two eV higher 
than the one for interior (bulk) bonds (5.2eV); similarly, 
edge bonds have a tendency to be shorter than bulk ones. 
As a consequence, the flake's interior undergoes a surface 
induced compression which is the more pronounced the 
smaller the flake is. This compression manifests itself 
in the way in which various observables depend on the 
flake size, N. For example, elastic constants (i.e. Lame 
parameters) of small flakes can exceed their bulk limit 
{fi + X ~ 70eV per ring, ^ « 51eV per ring, v — 0.162) 
by 30% {fi + A) or even by 70% (/x) . In comparison, the 
sheet (out of plane, buckling) stiffness, k « 1.2eV, is less 
sensitive to iV. Non-linearities remain weak (less than 
10% increase) as long as the ratio of out of plane ampli- 
tude and in plane wavelength of buckling is below 5%. 
To highlight the importance of quantum effects on elas- 
ticity we have also calculated the vibrational spectrum 
of graphene flakes. Quantum corrections affect mostly 
the sheet stiffness, k, lowering it significantly, about 26% 
within our DFT framework. 

Finally, based on these results we predict a pronounced 
shift of the Raman G- and D-peaks with decreasing 
fiake size to higher values. It is a natural consequence 
of the edge induced flake compression. The associated 
Griineisen parameters are 7d ~ 2.6eV/A^ and 7g = 
2.2eV/l2. 
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